###################################
# Analysis.R.
#
# Estimate treatment effect from 2024 Houston messaging campaign for
# Kousser, Thad, Seth J. Hill, Neil Malhotra, and Robert M. Stein. 2026. "Can a Public Information Campaign Increase Trust in American Elections?" Science Advances.
# Input file: Kousser_Hill_Malhotra_Stein_Data.csv
###################################


rm(list=ls())
if (.Platform$OS == "windows") { # Set working directory in location of script.
  .doit <- function() {
  frame_files <- lapply(sys.frames(), function(x) x$ofile);frame_files <- Filter(Negate(is.null), frame_files)
  PATH <- dirname(frame_files[[length(frame_files)]]); setwd(PATH) ; rm(PATH,frame_files)
  }
  try(.doit(),silent=T)
}

library(bit64)
library(data.table)
library(stargazer)
library(ordinal)  # for clm() function
library(sandwich)
library(lmtest)
options(digits=2)

save_and_display_stargazer <- function(models, file_name, keep_stats = c("n"), ...) {
  # Display on screen
  stargazer(models, ..., star.cutoffs = NA, type = "text", keep.stat = keep_stats)  
  # Write to HTML but not to screen.
  capture.output(
    stargazer(models, type = "html", ..., star.cutoffs = NA, out = file_name, keep.stat = keep_stats), file=NULL
  ) 
  flush.console()
}


###################
# Call in data.
###################
DT = fread("Kousser_Hill_Malhotra_Stein_Data.csv")

###################
# Linear probability moodels.
###################

# Functions to exectue regressions and calculate pre-registered standard errors.
# Estimate regressions and cluster SEs.
doLms <- function(DT, predictors=c("treated")) {
  # Models for hypotheses H1-H4 (Linear Probability Models)
  # H1: Confidence in personal vote count
  form = reformulate(termlabels=predictors,response="numeric_personal_conf")
  model_h1 <- lm(form, data = DT)

  # H2: Confidence in Harris County vote count
  form = reformulate(termlabels=predictors,response="numeric_harris_conf")
  model_h2 <- lm(form, data = DT)

  # H3: Trust in vote-by-mail
  form = reformulate(termlabels=predictors,response="numeric_vbm_trust")
  model_h3 <- lm(form, data = DT)

  # H4: Perception of vote fraud frequency
  form = reformulate(termlabels=predictors,response="numeric_fraud_freq")
  model_h4 <- lm(form, data = DT)

  # H5: Turnout in November 2024 election (controls for past turnout)
  form = reformulate(termlabels=c(predictors,"voted_22","voted_20"),response="voted_24")
  model_h5 <- lm(form, data = DT)

  # Compute clustered standard errors.
  library(sandwich)
  library(lmtest)
  clustered_ses <- list(
    sqrt(diag(vcovCL(model_h1, cluster = ~ household_id))),
    sqrt(diag(vcovCL(model_h2, cluster = ~ household_id))),
    sqrt(diag(vcovCL(model_h3, cluster = ~ household_id))),
    sqrt(diag(vcovCL(model_h4, cluster = ~ household_id)))
  )

  # Calculate MDEs for each model (using one-sided tests)
  mdes <- sapply(list(model_h1, model_h2, model_h3, model_h4), 
                 function(x) calculate_mde(x, one_sided = TRUE))

  # Output is a list of lists.
  return(list(models=list(model_h1, model_h2, model_h3, model_h4), mdes=mdes, turnout_model=model_h5))
}

# Minimum detectable effects function.
calculate_mde <- function(model, alpha = 0.05, power = 0.8, one_sided = TRUE) {
  # Get sample size
  n <- nobs(model)    
  # Get proportion treated
  p <- mean(model$model$treated)    
  # Get residual standard deviation
  sigma <- sigma(model)    
  # Critical values - adjust based on one_sided parameter
  if (one_sided) {
    t_alpha <- qt(1 - alpha, df = n - 2)
  } else {
    t_alpha <- qt(1 - alpha/2, df = n - 2)
  }
  t_beta <- qt(power, df = n - 2)    
  # Calculate and return MDE
  mde <- (t_alpha + t_beta) * sigma * sqrt(1/(n*p*(1-p)))    
  return(mde)
}

# Pre-registereed linear hypothesis
res <- doLms(DT)

  # Save out.
  save_and_display_stargazer(res$models, title = "Linear Probability Models", covariate.labels = c("Treatment"), dep.var.labels = c("Personal Vote Conf", "Harris Co Vote Conf", "Vote-by-Mail Trust", "Frequency Fraud"), add.lines = list(c("MDE (80% power, one-sided)", sprintf("%.3f", res$mdes))), notes = "Standard errors clustered on household",file_name = "Table-LPMs.html")
  # Separate table for turnout.
  save_and_display_stargazer(res$turnout_model, se = list(sqrt(diag(vcovCL(res$turnout_model, cluster = ~ household_id)))), title = "Turnout 2024", covariate.labels = c("Treatment","Turnout 2022","Turnout 2020"), dep.var.labels = c("Turnout"), add.lines = list(c("MDE (80% power, one-sided)", sprintf("%.3f", calculate_mde(res$turnout_model, one_sided = TRUE)))), notes = "Standard errors clustered on household", file_name = "Table-Turnout.html")


###################
# Reviewer-requested (not preregistered) hypothesis that treatment effects might
# interact with pre-treatment trust. Use machine learning estimated trust as
# moderator of treatment.
###################
# Break imputed non-citizen fraud into terciles.
cuts = DT[,quantile(imputed_noncitizen_fraud,probs=c(.33,.67))]
DT[,bottom_tercile_fraud := 1*(imputed_noncitizen_fraud <= cuts[1])]
DT[,top_tercile_fraud := 1*(imputed_noncitizen_fraud > cuts[2])]
res_zaller <- doLms(DT,predictors=c("treated*bottom_tercile_fraud","treated*top_tercile_fraud"))

  # Save out.
  save_and_display_stargazer(res_zaller$models, title = "Linear Probability Models", covariate.labels = c("Treatment","Bottom Third Imputed Fraud Frequency","Top Third Imputed Fraud Frequency","Treatment*Bottom Third","Treatment*Top Third"), dep.var.labels = c("Personal Vote Conf", "Harris Co Vote Conf", "Vote-by-Mail Trust", "Frequency Fraud"), add.lines = list(c("MDE (80% power, one-sided)", sprintf("%.3f", res_zaller$mdes))), notes = "Standard errors clustered on household",file_name = "Table-LPMs-Zaller.html")
  # Separate table for turnout.
  save_and_display_stargazer(res_zaller$turnout_model, se = list(sqrt(diag(vcovCL(res_zaller$turnout_model, cluster = ~ household_id)))), title = "Turnout 2024", covariate.labels = c("Treatment","Bottom Third Imputed Fraud Frequency","Top Third Imputed Fraud Frequency","Turnout 2022","Turnout 2020","Treatment*Bottom Third","Treatment*Top Third"), dep.var.labels = c("Turnout"), add.lines = list(c("MDE (80% power, one-sided)", sprintf("%.3f", calculate_mde(res_zaller$turnout_model, one_sided = TRUE)))), notes = "Standard errors clustered on household", file_name = "Table-Turnout-Zaller.html")

###################
# Machine learning imputation validation check
# (not pre-registered).
###################
# Predict each survey outcome with imputed fraud variable.
res_validation <- doLms(DT,predictors=c("imputed_noncitizen_fraud"))

  # Save out.
  save_and_display_stargazer(res_validation$models, title = "Linear Probability Models", covariate.labels = c("Machine learning imputed fraud belief"), dep.var.labels = c("Personal Vote Conf", "Harris Co Vote Conf", "Vote-by-Mail Trust", "Frequency Fraud"), notes = "Standard errors clustered on household", file_name = "Table-LPMs-Validation.html")

###################
# Create ordered factors for ordered logit.
###################
  DT[, fac_personal_conf := factor(y_personal_conf,levels=c("Not at all confident", "Not too confident", "Unsure", "Somewhat confident", "Very confident"))]
  DT[, fac_harris_conf := factor(y_harris_conf,levels=c("Not at all confident", "Not too confident", "Unsure", "Somewhat confident", "Very confident"))]
  DT[, fac_vbm_trust := factor(y_vbm_trust, levels=c("Don’t trust at all", "Trust a little", "Unsure", "Trust somewhat", "Trust a lot"))]
  DT[, fac_fraud_freq := factor(y_fraud_freq,levels=c("Not frequently at all","Somewhat infrequently", "Unsure", "Somewhat frequently", "Very frequently"))]

###################
# Ordered logit models mirroring H1-H4
###################
  # H1: Confidence in personal vote count
  ologit_h1 <- clm(fac_personal_conf ~ treated, data = DT)
  # H2: Confidence in Harris County vote count
  ologit_h2 <- clm(fac_harris_conf ~ treated, data = DT)
  # H3: Trust in vote-by-mail
  ologit_h3 <- clm(fac_vbm_trust ~ treated, data = DT)
  # H4: Perception of vote fraud frequency
  ologit_h4 <- clm(fac_fraud_freq ~ treated, data = DT)
  # Create stargazer table for ordered logit results
  ologit_models <- list(ologit_h1, ologit_h2, ologit_h3, ologit_h4)

  # Compute clustered standard errors.
  clustered_ol_ses <- list(
    sqrt(diag(vcovCL(ologit_h1, cluster = ~ household_id))),
    sqrt(diag(vcovCL(ologit_h2, cluster = ~ household_id))),
    sqrt(diag(vcovCL(ologit_h3, cluster = ~ household_id))),
    sqrt(diag(vcovCL(ologit_h4, cluster = ~ household_id)))
  )

  # Save out.
  save_and_display_stargazer(ologit_models, se = clustered_ol_ses, title = "Ordered Logit Models", covariate.labels = c("Treatment"), dep.var.labels = c("Personal Vote Conf", "Harris Co Vote Conf", "Vote-by-Mail Trust", "Frequency Fraud"), notes = "Standard errors clustered on household", model.names = FALSE, digits = 3, align = TRUE, file_name = "Table-Ordered-Logit.html")


###################
# Attrition check.
###################
# Questions were asked in a fixed order, so we have declining
# response rates to each question. Create variables noting
# attrition on each qx.
qxs = c("numeric_vbm_trust", "numeric_fraud_freq", "numeric_personal_conf", "numeric_harris_conf")
new_cols = gsub("numeric_","answered_",qxs)
DT[,(new_cols) := lapply(.SD,function(x) 1*(!is.na(x))), .SDcols = qxs]
print(t(DT[,lapply(.SD,sum),.SDcols=new_cols]))

  # OLS predicting whether respondents answered each question.
  attr_01 <- lm(answered_vbm_trust ~ treated, data=DT)
  attr_02 <- lm(answered_fraud_freq ~ treated, data=DT)
  attr_03 <- lm(answered_personal_conf ~ treated, data=DT)
  attr_04 <- lm(answered_harris_conf ~ treated, data=DT)

  # Compute clustered standard errors.
  clustered_attr_ses <- list(
    sqrt(diag(vcovCL(attr_01, cluster = ~ household_id))),
    sqrt(diag(vcovCL(attr_02, cluster = ~ household_id))),
    sqrt(diag(vcovCL(attr_03, cluster = ~ household_id))),
    sqrt(diag(vcovCL(attr_04, cluster = ~ household_id)))
  )

  # Save out.
  save_and_display_stargazer(attr_01, attr_02, attr_03, attr_04, se = clustered_attr_ses, title = "Attrition", covariate.labels = c("Treatment"), dep.var.labels = c("Vote By Mail Trust", "Vote Fraud Frequency", "Personal Vote Confidence","Harris Vote Confidence"), notes = "Standard errors clustered on household", model.names = FALSE, digits = 3, align = TRUE, file_name = "Table-Attrition.html")



###################
# Turnout rates reported in paper.
###################
cat("\n2024 Turnout rates:\nFull target population:\n")
  print(DT[,mean(voted_24)])
cat("For those who answered the first post-treatment survey question on vote by mail:\n")
  print(DT[answered_vbm_trust == 1,mean(voted_24)])


###################
# Post hoc: Voting mode.
###################
# Condition on having turned out.
DT[,mail_24 := ifelse(turnout_mode_24g == "MAIL-IN",1,ifelse(turnout_mode_24g == "",NA,0))]
DT[,table(turnout_mode_24g,mail_24,exclude=NULL)]
DT[,mail_22 := 1*(vh_22_g == "Voted Absentee/by Mail")]
DT[,mail_20 := 1*(vh_20_g == "Voted Absentee/by Mail")]

  mode_lm <- lm(mail_24 ~ treated, data = DT)
  # Mail 2022 or 2020 voters.
  mode2_lm <- lm(mail_24 ~ treated, data = DT[mail_22 == 1 | mail_20 == 1,])
  # In-person or early voters.
  mode3_lm <- lm(mail_24 ~ treated, data = DT[mail_22 == 0 & mail_20 == 0,])

  # Compute clustered standard errors.
  clustered_md_ses <- list(
     sqrt(diag(vcovCL( mode_lm, cluster = ~ household_id)))
    ,sqrt(diag(vcovCL(mode2_lm, cluster = ~ household_id)))
    ,sqrt(diag(vcovCL(mode3_lm, cluster = ~ household_id)))
  )

  # Save out.
  save_and_display_stargazer(mode_lm, mode2_lm, mode3_lm, se=clustered_md_ses, title = "Turnout Mode 2024", covariate.labels = c("Treatment"), column.labels = c("All 2024 Voters","By Mail 2020 or 2022","Not By Mail"), dep.var.labels = c("Voted by Mail 2024","","(Voted 2024)"), notes = "Standard errors clustered on household", model.names = FALSE, digits = 3, align = TRUE, file_name = "Table-Voting-Mode.html")


###################
# Post-Hoc: Balance table for treatment vs control for each
# survey question in anlaysis.
###################

# Function to check balance b/w two data.tables.
check_balance <- function(treatment_dt, control_dt) {
  # Ensure both data.tables have the same columns
  if (!all(names(treatment_dt) == names(control_dt))) {
    stop("Treatment and control data.tables must have the same columns")
  }
  
  # Initialize result data.table
  result <- data.table(Variable = character(), Group = character(), Summary = character())
  
  for (var in names(treatment_dt)) {
    if (is.numeric(treatment_dt[[var]])) {
      # For numeric variables, calculate mean and standard deviation
      treat_mean <- mean(treatment_dt[[var]], na.rm = TRUE)
      treat_sd <- sd(treatment_dt[[var]], na.rm = TRUE)
      control_mean <- mean(control_dt[[var]], na.rm = TRUE)
      control_sd <- sd(control_dt[[var]], na.rm = TRUE)
      
      result <- rbind(result,
                      data.table(Variable = var, Group = "Treatment", Summary = sprintf("%.2f (%.2f)", treat_mean, treat_sd)),
                      data.table(Variable = var, Group = "Control", Summary = sprintf("%.2f (%.2f)", control_mean, control_sd)))
    } else {
      # For categorical variables, calculate proportions
      treat_table <- table(treatment_dt[[var]])
      control_table <- table(control_dt[[var]])
      
      # Function to replace blank categories and calculate proportions
      replace_blank_and_prop <- function(tbl) {
        names(tbl) <- ifelse(names(tbl) %in% c("", " "), "(BLANK)", names(tbl))
        prop <- prop.table(tbl)
        paste(names(prop), sprintf("%.1f%%", prop * 100), collapse = ", ")
      }
      
      result <- rbind(result,
                      data.table(Variable = var, Group = "Treatment", Summary = replace_blank_and_prop(treat_table)),
                      data.table(Variable = var, Group = "Control", Summary = replace_blank_and_prop(control_table)))
    }
  }
  
  # Set key for easy viewing
  setkey(result, Variable, Group)
  
  return(result)
}

## Covariates on which to check balance.
DT[,Female := 1*(sex == "F")]
covars = c("age", "Female", "voted_24", "voted_22", "voted_20")

# Create balance table.
  # H1: Confidence in personal vote count
  BT = check_balance(DT[treated == 1 & answered_personal_conf == 1,covars,with=F], DT[treated == 0 & answered_personal_conf == 1,covars,with=F])
  # Make wide.
  WD = dcast(BT,Variable ~ Group,value.var="Summary")
  WD[,Question := "Personal Vote Conf"]

  # H2: Confidence in Harris County vote count
  BT = check_balance(DT[treated == 1 & answered_harris_conf == 1,covars,with=F], DT[treated == 0 & answered_harris_conf == 1,covars,with=F])
  WD1 = dcast(BT,Variable ~ Group,value.var="Summary")
  WD1[,Question := "Harris Co Vote Conf"]
  WD = rbind(WD,WD1,use.names=T)

  # H3: Trust in vote-by-mail
  BT = check_balance(DT[treated == 1 & answered_vbm_trust == 1,covars,with=F], DT[treated == 0 & answered_vbm_trust == 1,covars,with=F])
  WD1 = dcast(BT,Variable ~ Group,value.var="Summary")
  WD1[,Question := "Vote-by-Mail Trust"]
  WD = rbind(WD,WD1,use.names=T)

  # H4: Perception of vote fraud frequency
  BT = check_balance(DT[treated == 1 & answered_fraud_freq == 1,covars,with=F], DT[treated == 0 & answered_fraud_freq == 1,covars,with=F])
  WD1 = dcast(BT,Variable ~ Group,value.var="Summary")
  WD1[,Question := "Frequency Fraud"]
  WD = rbind(WD,WD1,use.names=T)

  # Also include, for reference, full target population numbers.
  BT = check_balance(DT[treated == 1,covars,with=F], DT[treated == 0,covars,with=F])
  WD1 = dcast(BT,Variable ~ Group,value.var="Summary")
  WD1[,Question := "Full Target Population"]
  WD = rbind(WD,WD1,use.names=T)

  # Rename variables.
  WD[Variable == "age",Variable := "Age"]
  WD[Variable == "voted_20",Variable := "Turnout 2020"]
  WD[Variable == "voted_22",Variable := "Turnout 2022"]
  WD[Variable == "voted_24",Variable := "Turnout 2024"]
  setcolorder(WD,c("Variable","Question"))
  setkey(WD,Variable,Question)

  print(WD)
  fwrite(WD,"BalanceBySurveyQuestion.csv")


###################
# Differential attrition observed. Calculate Lee bounds following
# pre-analysis plan.
###################

lee_bounds <- function(data, treatment_var, outcome_var, attrition_var, n_bootstrap = NULL, conf_level = 0.95, verbose = TRUE) {
  # Helper function to calculate bounds
  calc_bounds <- function(data) {
    # Split observed data into treatment and control groups
    observed <- data[data[[attrition_var]] == 1, ]
    t_df <- observed[observed[[treatment_var]] == 1, ]
    c_df <- observed[observed[[treatment_var]] == 0, ]
    
    # Calculate attrition rates
    attrition_rate_treated <- sum(data[[treatment_var]] == 1 & data[[attrition_var]] == 0) / sum(data[[treatment_var]] == 1)
    attrition_rate_control <- sum(data[[treatment_var]] == 0 & data[[attrition_var]] == 0) / sum(data[[treatment_var]] == 0)
    
    # Determine trimming rate
    trim_rate <- abs(attrition_rate_treated - attrition_rate_control)
    
    # Sort by outcome to trim
    treated_sorted <- t_df[order(t_df[[outcome_var]]), ]
    control_sorted <- c_df[order(c_df[[outcome_var]]), ]

    # Trim for differential attrition.
    if (attrition_rate_treated > attrition_rate_control) {
      # Trim control group to match observed attrition rate in treated group.
      n_trim <- floor(nrow(control_sorted) * trim_rate)
      if (n_trim >= nrow(control_sorted)) return(c(NA, NA))
      treated_trimmed <- treated_sorted
      control_trimmed <- control_sorted[-seq_len(n_trim), ]
    } else {
      # Trim treatment group to match observed attrition rate in control group.
      n_trim <- floor(nrow(treated_sorted) * trim_rate)
      if (n_trim >= nrow(treated_sorted)) return(c(NA, NA))
      treated_trimmed <- treated_sorted[-seq_len(n_trim), ]
      control_trimmed <- control_sorted
    }
    
    # Calculate bounds
    lower_bound <- mean(treated_trimmed[[outcome_var]], na.rm = TRUE) - mean(control_trimmed[[outcome_var]], na.rm = TRUE)
    upper_bound <- mean(treated_sorted[[outcome_var]], na.rm = TRUE) - mean(control_sorted[[outcome_var]], na.rm = TRUE)
    c(lower_bound, upper_bound)
  }
  
  # Calculate initial bounds
  initial_bounds <- calc_bounds(data)
  
  # Bootstrap for inference
  if (!is.null(n_bootstrap)) {
    bootstrap_bounds <- replicate(n_bootstrap, {
      bootstrap_data <- data[sample(seq_len(nrow(data)), replace = TRUE), ]
      calc_bounds(bootstrap_data)
    })
  
    # Remove invalid (NA) results
    bootstrap_bounds <- bootstrap_bounds[, complete.cases(bootstrap_bounds)]
    
    # Confidence intervals
    alpha <- 1 - conf_level
    lower_ci <- apply(bootstrap_bounds, 1, function(x) quantile(x, probs = alpha / 2, na.rm = TRUE))
    upper_ci <- apply(bootstrap_bounds, 1, function(x) quantile(x, probs = 1 - alpha / 2, na.rm = TRUE))
  } else {
    lower_ci <- upper_ci <- NA
  }
  
  # Print results if verbose
  if (verbose) {
    cat("Results of Lee Bounds for outcome",outcome_var,":\n")
    cat("Lower Bound:", initial_bounds[1], "\n")
    cat("Upper Bound:", initial_bounds[2], "\n")
    if (!is.null(n_bootstrap)) {
      cat(conf_level * 100, "% CI for Lower Bound:", lower_ci[1], "to", upper_ci[1], "\n")
      cat(conf_level * 100, "% CI for Upper Bound:", lower_ci[2], "to", upper_ci[2], "\n")
    }
  }
  
  # Return results as a list
  list(
    bounds = initial_bounds,
    ci_lower = lower_ci,
    ci_upper = upper_ci
  )
}

# Bootstrap confidence intervals (set to NULL for none).
n_boots = NULL
n_boots = 50

# Models for hypotheses H1-H4 (Linear Probability Models)
set.seed(2025)
  # H1: Confidence in personal vote count
  lee_h1 <- lee_bounds(data = DT, treatment_var="treated", attrition_var = "answered_personal_conf", n_bootstrap = n_boots, outcome_var = "numeric_personal_conf")

  # H2: Confidence in Harris County vote count
  lee_h2 <- lee_bounds(data = DT, treatment_var="treated", attrition_var = "answered_harris_conf", n_bootstrap = n_boots, outcome_var = "numeric_harris_conf")

  # H3: Trust in vote-by-mail
  lee_h3 <- lee_bounds(data = DT, treatment_var="treated", attrition_var = "answered_vbm_trust", n_bootstrap = n_boots, outcome_var = "numeric_vbm_trust")

  # H4: Perception of vote fraud frequency
  lee_h4 <- lee_bounds(data = DT, treatment_var="treated", attrition_var = "answered_fraud_freq", n_bootstrap = n_boots, outcome_var = "numeric_fraud_freq")

  # H5: Turnout in November 2024 election (controls for past turnout).
  # Use answered first question on survey (vbm) as attrition measure.
  lee_h5 <- lee_bounds(data = DT, treatment_var="treated", attrition_var = "answered_vbm_trust", n_bootstrap = n_boots, outcome_var = "voted_24")

# Turn to list of results to pass to stargazer.
lee_bounds_table <- function(results_list, model_names = NULL, conf_level = 0.95, verbose = FALSE) {
  # Check if results_list is a list of lee_bounds outputs
  if (!all(sapply(results_list, function(res) is.list(res) && all(c("bounds", "ci_lower", "ci_upper") %in% names(res))))) {
    stop("results_list must be a list of results from the lee_bounds function.")
  }
  
  # Generate default model names if not provided
  if (is.null(model_names)) {
    model_names <- paste("Model", seq_along(results_list))
  }
  
  # Initialize a data frame to hold the results
  table_data <- data.frame(
    Model = model_names,
    Lower_Bound = sapply(results_list, function(res) res$bounds[1]),
    Upper_Bound = sapply(results_list, function(res) res$bounds[2]),
    CI_Lower_Bound = sapply(results_list, function(res) res$ci_lower[1]),
    CI_Upper_Bound = sapply(results_list, function(res) res$ci_upper[2])
  )
  
  # Verbose printing of the table
  if (verbose) {
    cat("\nLee Bounds Summary Table\n")
    cat(sprintf("Confidence Level: %.2f%%\n", conf_level * 100))
    cat("-----------------------------------------\n")
    cat(sprintf("%-10s | %10s | %10s | %10s | %10s\n", 
                "Model", "Lower Bound", "Upper Bound", "CI Lower", "CI Upper"))
    cat("-----------------------------------------\n")
    for (i in seq_len(nrow(table_data))) {
      cat(sprintf("%-10s | %10.3f | %10.3f | %10.3f | %10.3f\n",
                  table_data$Model[i], 
                  table_data$Lower_Bound[i],
                  table_data$Upper_Bound[i],
                  table_data$CI_Lower_Bound[i],
                  table_data$CI_Upper_Bound[i]))
    }
    cat("-----------------------------------------\n")
  }
  
  # Return the table data for use with stargazer
  return(table_data)
}

  results_table <- lee_bounds_table(list(lee_h1, lee_h2, lee_h3, lee_h4, lee_h5), model_names = c("Personal Vote Conf", "Harris Co Vote Conf", "Vote-by-Mail Trust", "Frequency Fraud", "Turnout"))

  # Save out.
  save_and_display_stargazer(results_table, summary = FALSE, title = "Lee Bounds Summary Table", digits = 3, file_name = "Table-Lee-Bounds.html")
